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U ' We present and analyse in this paper three models of coupled continuum equations all united by 

I ' a common theme: the intuitive notion that sandpile surfaces are left smoother by the propagation 

'■^ . of avalanches across them. Two of these concern smoothing at the 'bare' interface, appropriate to 

(Ih ' intermittent avalanche flow, while one of them models smoothing at the effective surface defined 

O I by a cloud of flowing grains across the 'bare' interface, which is appropriate to the regime where 

I ^1 ■ avalanches flow continuously across the sandpile. 

,__! ■ PACS NOS.: 05.40+j, 05.70.Ln, 46.10.+Z, 64.60.Ht 

cn ■ 

oo ; 

O . I. INTRODUCTION 

(N; 

r ■ The dynamics of sandpiles have intrigued researchers in physics over recent years [0,0] with a great deal of effort 

^s^ being devoted to the development of techniques involving for instance cellular automata p,H , continuum equations 
*>••; P^nl a-iid Monte Carlo schemes Isl to investigate this very complex subject. However what have often been lost sight 

C^ of in all this complexity are some of the extremely simple phenomena that are exhibited by granular media which still 
G , remain unexplained. 

I ' One such phenomenon is that of the smoothing of a sandpile surface after the propagation of an avalanche [BJ. 

' ^ It is clear what happens physically: an avalanche provides a means of shaving off roughness from the surface of a 

^ sandpile by transferring grains from bumps to available voids ||Q|, and thus leaves in its wake a smoother surface. 

fx \ However, surprisingly, researchers have not to our knowledge come up with models of sandpiles that have exhibited 

• • . this behaviour. 

In particular what has not attracted enough attention in the literature is the qualitative difference between the 

situations which obtain when sandpiles exhibit intermittent and continuous avalanches pO| . In this paper we examine 

^ ■ both the latter situations, via distinct models of sandpile surfaces. 

5^ , A particular experimental paradigm that we choose to put our discussions in context is that of sand in rotating 

cylinders |ll| , |l2(| . In the case when sand is rotated slowly in a cylinder, intermittent avalanching is observed; thus 

sand accumulates in part of the cylinder to beyond its angle of repose [Esl and is then released via an avalanche 
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process across the slope. This happens intermittently, since the rotation speed is less than the characteristic time 
between avalanches. By contrast, when the rotation speed exceeds the time between avalanches, we see continuous 
avalanching on the sandpile surface. Though this phenomenon has been observed ||l^ and analysed physically [lC| ] 
in terms of avalanche statistics, we are not aware of measurements which measure the characteristics of the resulting 
surface in terms of its smoothness or otherwise. 

What we focus on here is precisely this aspect, and make predictions which we hope will be tested experimentally. 
In order to discuss this, we introduce first the notion that granular dynamics is well described by the competition 
between the dynamics of grains moving independently of each other and that of their collective motion within clusters 
H . A convenient way of representing this is via coupled continuum equations with a specific coupling between mobile 
grains p and clusters h on the surface of a sandpile g . In the regime of intermittent avalanching, we expect that the 
interface will be the one defined by the 'bare' surface, i.e. the one defined by the relatively immobile clusters across 
which grains flow intermittently. This then implies that the roughening characteristics of the h profile should be 
examined. The simplest of the three models we discuss in this paper (an exactly solvable model referred to hereafter 
as Case A) as well as the most complex one (referred to hereafter as Case C) treat this situation, where we obtain 
in both cases an asymptotic smoothing behaviour in h. When on the other hand, there is continuous avalanching, 
the flowing grains provide an effective film across the bare surface and it is therefore the species p which should be 
analysed for spatial and temporal roughening. In the model hereafter referred to as Case B we look at this situation, 
and obtain the surprising result of a gradual crossover between purely diffusive behaviour and hypersmooth behaviour. 
In each case we present analytical results pertaining to the continuum models and compare the predictions so obtained 
with the results obtained by numerical simulations of the corresponding discretised equations. 

In general, the complexity of sandpile dynamics leads us to equations which are coupled, nonlinear and noisy: these 
equations present challenges to the theoretical physicist in more ways than the obvious ones to do with their detailed 
analysis and/or their numerical solutions. In particular, our analysis of Case C reveals the presence of hidden length 
scales whose existence was suspected analytically, but not demonstrated numerically in earlier work P,M . 

The normal procedure for probing temporal and spatial roughening in interface problems is to determine the 
asymptotic behaviour of the interfacial width with respect to time and space, via the single Fourier transform. Here 
only one of the variables, {x, t) is integrated over in Fourier space, and appropriate scaling relations are invoked to 
determine the critical exponents which govern this behaviour. However, it turns out that this leads to ambiguities for 
those classes of problems where there is an absence of simple scaling, or to be more specific, where multiple length 
scales exist. In such cases we demonstrate that the double Fourier transform (where both time and space are integrated 
over) yields insights that are harder to obtain via the single Fourier transform. 

This point is illustrated by Case A, an exactly solvable model that we introduce; we then use it to understand Case 
C, a nonlinear model where our analytical results are clearly only approximations to the truth. 

In order to make some of these ideas more concrete, we now review some general facts about rough interfaces [[l5| . 
Three critical exponents, a, /3, and z, characterise the spatial and temporal scaling behaviour of a rough interface. 
They are conveniently defined by considering the (connected) two-point correlation function of the heights 

S{x -x',t~ t') = {h{x, t)h{x', t')) - {h{x, t)){h{x', t')). (1) 

We have 

S'(x,0) - |a;p" (l^l^oo) and S'(0,<) - jip'^ (|t| ^ c5o), 

and more generally 

S{x,t)^\x\'"F{\t\/\x\^) 

in the whole long-distance scaling regime (x and t large). The scaling function F is universal in the usual sense; a 
and z — a/P are respectively referred to as the roughness exponent and the dynamical exponent of the problem. In 
addition, we have for the full structure factor which is the double Fourier transform S{k,uj) 

which gives in the limit of small k and ui, 

S{k, w = 0) - k-'^-^"-' {k -^ 0) and S(k = 0, cj) - w^i-^/J-i/z ^^ ^ q^ ^2) 

The scaling relations for the corresponding single Fourier transforms are 

S'(fc,i = 0)~fc~i-2a (;,_o) and S {x = 0,lj) ^ uj-'^-^^ {uj ^ 0) (3) 



In particular we note that the scahng relations for S{k,u!) (Eq.(g)) always involve the simultaneous presence of a 
and P, whereas those corresponding to S(x,uj) and S(k,t) involve these exponents individually. Thus, in order to 
evaluate the double Fourier transforms, we need in each case information from the growing as well as the saturated 
interface (the former being necessary for /3 and the latter for a) whereas for the single Fourier transforms, we need only 
information from the saturated interface for S{k, t — 0) and information from the growing interface for S{x = 0,lu). On 
the other hand, the information that we will get out of the double Fourier transform will provide a more unambiguous 
picture in the case where multiple length scales are present, something which cannot easily be obtained in every case 
with the single Fourier transform. 

In Sections II, III and IV we present, analyse and discuss the results of Cases A, B and C respectively. Finally, 
in Section V, we reflect on the unifying features of these models, and make some educated guesses on the dynamical 
behaviour of real sandpile surfaces. 

II. CASE A: THE EDWARDS- WILKINSON EQUATION WITH FLOW 

Our first model involves a pair of linear coupled equations, where the equation governing the evolution of clusters 
("stuck" grains) h is closely related to the very well-known Edwards- Wilkinson (EW) model [|l6[. The equations are: 

dh(x t) 

^ ' ' ^ Dh\7'^h{x,t) + cS/h{x,t)+r]{x,t) (4a) 



dt 

DpV^p{x, t) ~ cVh{x, t) (4b) 



dp{x, t) _ ^2 
dt -^" 



where the first of the equations describes the height h{x, t) of the sandpile surface at (x, t) measured from some mean 
(/i), and is precisely the EW equation in the presence of the flow term cV/i. The second equation describes the 
evolution of flowing grains, where p{x,t) is the local density of such grains at any point {x,t). As usual, the noise 
r]{x,t) is taken to be Gaussian so that: 

{rj{x, t)rj{x',t')) = A^S{x - x')5{t - t'). 

with A the strength of the noise. Here, (• • •) refers to an average over space as well as over noise. 

A. Analysis of the decoupled equation in h 

For the purposes of analysis, we focus on the first of the two coupled equations (Eq.(Ea)) presented above, 

dh 

— ^ DhV^h + c\/h + T]{x,t) 

noting that this equation is essentially decoupled from the second. (This statement is, however, not true in reverse, 
which has implications to be discussed later). We note that this is entirely equivalent to the Edwards- Wilkinson 
equation |16| in a frame moving with velocity c 



x' = X + ct , t' = t 

and would on these grounds expect to find only the well-known EW exponents a = 0.5 and j3 — 0.25 |l|]. This would 
be verified by naive single Fourier transform analysis of Eq. ( [4a|) which yields these exponents via Eq. (| 
Equation ( pa| ) can be solved exactly as follows. The propagator G{k,Lu) is 

Gh{k,uj) = {~iuj + Dhk'^ + ikcy^ 

This can be used to evaluate the structure factor 

{h{k,uj)h{k\uj')) 



Sh{k,uj) = 



6{k + k')SiLu + uj') 



which is the Fourier transform of the full correlation function Sh{x — x',t — t') defined by Eq.(n|). The solution for 
Sh{k,u}) so obtained is: 



This is illustrated in Fig.|l| while representative graphs for Sh{k,uj = 0) and Sh{k = 0,w) are presented in Figs.g 
and respectively. Before proceeding further, we make the following observation about the double Fourier transform 
5/1 (fc = 0, w); this shows an w~^ behaviour coming from Eq.(^), which we will also see later. We mention here that the 
ubiquity of this uj~^ arises from the form of the scaling relation Eq.(0), which is relevant for frequencies uj < lJc ~ k^'^ , 



whereas for to > uJc the high frequency behaviour takes over giving cj^^ (cf. Eq.(A2) in the Appendix). As fc = for 
the purposes of calculation of this structure factor, it is always the high frequency behaviour that dominates, leading 
to the ubiquity of oj"^ whenever it is measured. 

It is obvious from Eq.(P) that Sh{k,uj) does not show simple scaling. More explicitly, if we write 

with fco = c/Dfi, and luq — c^ /Dh, we see that there are two limiting cases : 

for / 
Eqs 



k^' 



• for k 3> fco, Sf^ ^{k, cj = 0) ~ k^; using again Sf^^{k = 0, uj) ^ u)^, we obtain a^ — 1/2 and f3h — 1/4, z^ — 2 via 



• for k <^ fco, 5,j ^(fc, UJ — 0) ^ k^; using the fact that the limit Sf^^{k — 0, uj) is always w^, this is consistent with 
the set of exponents a^ — 0, f3h — and z^ = I via Eqs.(0). 

The first of these contains no surprises, being the normal EW fixed point |16| , while the second represents a new, 
'smoothing' fixed point. 

We now explain this smoothing fixed point via a simple physical picture. The competition between the two terms 
in Eq.(Bq) determines the nature of the fixed point observed: when the diffusive term dominates the flow term, the 
canonical EW fixed point is obtained, in the limit of large wavevectors fc. On the contrary, when the flow term 
predominates, the effect of diffusion is suppressed by that of a travelling wave whose net result is to penalise large 
slopes; this leads to the smoothing fixed point obtained in the case of small wavevectors k. We emphasise however, 
that this is a toy model of smoothing, which will be used to illuminate the discussion of models B and C below. 

B. Coupled equations: a model of smoothing 

We realise from the above that the interface h is smoothed because of the action of the flow term which penalises the 



sustenance of finite gradients \7h in Eq.([4aD. However, Eq.(|4a|) is effectively decoupled from Eq.(4b), while Eq.([4b|) is 
manifestly coupled to Eq.(Q). In order for the coupled Eqs. (4) to qualify as a valid model of sandpile dynamics, we 
would need to ensure that no instabilities are generated in either of these by the coupling term cV/i. 

In this spirit, we look first at the value of p averaged over the sandpile, as a function of time (Fig.pk). We observe 
that the incursions of {p) into negative values are limited to relatively small values, suggesting that the addition of a 
constant background of p exceeding this negative value would render the coupled system meaningful, at least to a first 
approximation. In order to ensure that this average does not involve wild fluctuations, we examine the fluctuations 
in yO, viz.-\/(/3^) — (p)^ (Fig.pp). The trends in that figure indicate that this quantity appears to saturate, at least 
upto computationally accessible times. Finally we look at the minimum and maxim,um value of p at any point in the 
pile over a large range of times (Fig.gp); this appears to be bounded by a modest (negative) value of 'bare' p. Our 
conclusions are thus that the fluctuations in p saturate at computationally accessible times and that the negativity of 
the fluctuations in p can always be handled by starting with a constant po, a constant 'background' of flowing grains, 
which is more positive than the largest negative fluctuation. 

Physically, then, the above implies that at least in the presence of a constant large density po of flowing grains, 
it is possible to induce the level of smoothing corresponding to the fixed point a = j3 = 0. This model is thus one 
of the simplest possible ways in which one can obtain a representation of the smoothing of the 'bare surface' that is 
frequently observed in experiments on real sandpiles after intermittent avalanche propagation 0. 



III. CASE B: A SIMPLE FORM OF COUPLING, WITH COMPLEX CONSEQUENCES 

Our model equations, first presented in |g| involve a simple coupling between the species h and p, where the 
transfer between the species occurs only in the presence of the flowing grains and is therefore relevant to the regime 
of continuous avalanching when the duration of the avalanches is large compared to the time between them. The 
equations are: 

dhix,t) _ „ _2 



dt 
dp{x, t) 



= DhV'h{x, t) - T{h, p) + 77ft (x, t) (6a) 

= DpV^p{x,t) + T{h,p) + T^p{x,t) (6b) 



dt 
T{h,p) = -pp{Wh) (6c) 

where the terms rjh{x,t) and rip(x,t) represent Gaussian white noise as usual: 

{Tjhix, t)'nh{x', t')) = ^l5{x - x')6it - t') 
{T^pix, t)7^p{x', t')) = AlS{x - x')S{t - t') 

and the (• • •) stands for average over space as well as noise. 

A simple physical picture of the coupling or 'transfer' term T{h, p) between h and p is the following: flowing grains 
are added in proportion to their local density to regions of the interface which are at less than the critical slope, and 
vice versa, provided that the local density of flowing grains is always non-zero. This form of interaction becomes zero 
in the absence of a finite density of flowing grains p (when the equations become decoupled) and is thus the simplest 
form appropriate to the situation of continuous avalanching in sandpiles. We analyse in the following the profiles of 
h and p consequent on this form. 

It turns out that a singularity discovered by Edwards iQ three decades ago in the context of fluid turbulence is 
present in models with a particular form of the transfer term T ; the above is one example, while another example is 
the model due to Bouchaud et al. (BCRE) |R] where 

T = -z/Vft. - pp{Vh) 

and the noise is present only in the equation of motion for h. This singularity, the so-called infrared divergence, 
largely controls the dynamics and produces unexpected exponents. 

A. Theoretical analysis 

We carry out first the theoretical analysis of Eqs.(6). An examination of the above equations reveals the presence of 
two likely length scales in each, one associated with the diffusive motion, and the other with the so-called transfer term 
T{h, p), representing the coupling between the two species. In these circumstances, a renormalisation group analysis 
would clearly be inappropriate due to the breakdown of simple scaling. In recent years, however, a self-consistent 
mode coupling analysis used hitherto in dynamic critical phenomena |l9|] has been used to look at in particular the 
Kardar-Parisi-Zhang (KPZ) equation fl^JSQ] and we extend its use to the case of the coupled equations presented 
here. 

In this method we set up equations (to one-loop order) for the correlation functions and self-energies in terms of the 
full Green's functions, correlation functions and vertices using assumed scaling forms for each. The critical exponents 
a and /3 defined above are obtained from the self-consistent solutions of these equations using Dh — Dp. 

Focusing on the h variable to start with, we define the Green's functions and the correlation functions of the h and 
p variables 

/ Sh{k,Lu) \ 1 

''^ ''^^ \6r]{k',n)/ S{k + k')Siu; + n) 

/ Spik,u;) \ 1 

"^ '""^ \Srj{k',n)/ S{k + k')S{u; + n) 

{h{k,uj)h{k',n)) 



Sh{k,uj) = 
Sp{k,uj) = 



S{k + k')5{io -h n) 

{p{k,uj)p{k',n)) 

6{k + k')S{io + n) 



The analysis of these functions wiU be in terms of a weak scahng hypothesis which states 



G.(fc,^) = fc--A(^,^) 



Gpik.Lo) 






A strong scaUng would imply the existence of a single time scale i.e. Zh 



As we show below, this cannot be the 



case here. The absence of strong scaling implies that the roughness exponents ah and Up may become functions of k. 
We consider the full Green's function Gh{k,uj), which is given via the well-known Dyson equation [gl| , 

G-\k,u;) = Gr\k,u;) + j:h{k,u;) 
Here, the zeroth order Green's function is 

The scaling forms of the functions Gh{k,uj) and Sh{k,Lo) are given by, in the limit tj ^ 0, 

1 



' + fc2 + fc^h 



1 



Gh{k,u) 
Sh{k,uj) 



Similar scaling relations hold for the species p. 

To one- loop order, the self-energy S;i(fc) is given by (Fig.^b) 



1 



OJ^ 



T.h{k,Lu) = Z-*^ / ^ / 2^^''(^ -1,^ - ^)Sp{q, 



n)k{k-q) 



dq 

2^ 



dn 

2^ 



i{ui — f2) -f E/i(fc — g, cj — rj) 



k[k ~ q) 



q 



l+2ap 



2Ep(g,fi) 



r!2 + |Sp(g,rj)| 



(7a) 
(7b) 



where the second line follows from the first in the limit of small fl. We note that due to the presence of the term 
q-i-2ap^ the integral is dominated by the singularity in the integrand at g — > 0. This 'infrared divergence' which 
results from the divergence of the internal momenta q, is very different from the usual divergences encountered in 
critical phenomena where the latter occur for small wave numbers and are associated with long wavelength instabilities 
in the external momenta. In this case due to the infrared divergence in the above equation in the internal momenta 
g, the integral diverges for any value of the external momenta k, so long as ap > 0. 

We thus need either to evaluate the integral with a lower cut-off fco or to introduce a suitable regulator. We follow 
the first of these procedures for the above equation, and the second of the procedures to do with the corresponding 
quantity, Sp{k,ijj), for p. 

We then proceed to evaluate the self-energy at zero external frequency, i.e T,h{k,uj = 0) from Eq.dTa). As g ^ 
we can approximate Gh{k — q, —Q.) by 

G~^(fc, ~n) ^in + k^ + Y^hik, -Vi) 
^e + Yh{k,Q) 

where the second line follows from the fact that we are looking at the q ~ limit of the internal frequency Q, ^ q^*^ . 
As E/j(fc, 0) ~ fc^*", the small k behaviour of Gh{k) is dominated by Y^h{k) for z/j < 2 , i.e 

The integral in Eq. ( |7£i| ) becomes in the limit of zero external frequencies 

p.^e 



^h{k) = 



Yh{k) 



Using the scaling form for the single Fourier transform (Eq.(@)) we find 

E,(fc) = M^fc2l],(fc)-iC,|g^ 



We now have to evaluate the integral by cutting off the momentum integration at fcp ^ 1 , i.e. we follow the first of 
the procedures given above to handle the infrared divergence. This gives, after some simplification, 

4:TTap 
From the above equation with the scaling relation T,h{k) ~ k^'' we find, on equating powers of k, 

Zh ^ I 

We note here that the presence of the term pVh could in principle cause the vertex fi to renormalise, leading to a 
correction to Zh- In these circumstances, the expression for the self-energy S;i(fc,cj ~ 0) is given by 



E;.(fc, c. = 0) = /i2 y J y ^T,{k, q, k - q)Gh{k - q, -n)Spiq, n)k{k - q) 



(8) 



where we have introduced a three-point vertex function r3(fc, q.k — q) in Eq.(^q). Assuming that as g ^ 0, we can 
write the asymptotic form for the three-point vertex as. 



T:,{k,q,k-q)^k-- 



(9) 



we find 



z/i = 1 



In the event that numerical results suggest z^ 7^ 1 we will have to incorporate this new renormalised vertex into our 
calculations. 

Next we examine the correlation function for /i, Sh{k^oj)^ which to one-loop order is given by (FigJqa) 



Sh{k,u) = 



c^2 + |E^(fc,a;)|2 
1 

TJ^+Wk^ 



.2 + |5]^(fc,u;)|2 



l + f^'l^l^\k-q\'S,{k-q,cj-n)S,{q,n) 

2Ep{q,rt) 



2 I dq f dn |fc — gp 1 



2tT J 2TT\k- g|l+2"/x gl+2"p \n^ + |I]^(g, f7)|2 

2E^(fc — q,LO — fl) 



1 + m' 



dq \k — q 



1-2q. 






2^ 9'+'"" W + {J:p{q) + ^hik-q)f 



(10a) 

(10b) 
(10c) 



The frequency-dependent self-energy I]^(/c,a;) in the above is given by evaluating the integral over the internal 
frequency H. in Eq.(7b). This leads to 



S/i(fc,w) w ^2 



2 f dq k{k - q) 



A 



2tt gi+2"p -iuj + Y.p{q) + ^h{k-q) 

2 7,-2Qp 



2 -A K Kq 



Anap —iuj + Tok 



TIP 



'iuj + Tok 



(11a) 
(lib) 
(lie) 



where Fq = /ifcg '' \ ^^a ' ^^'^ ^^^ second line in the above follows from taking a g — > limit and introducing a cutoff 

wavevector /cq in the integral on the first line. Introducing this expression for S/j(fc,ct') in Eq.(lOc) and recognising 
that the divergence due to g^(i+2"p) dominates the integral we find 



p4jL4 



Cj2 + r2/c2 



1 + ^fc-2"P 



fc2 



Fofc 



47rap'° /ci+2"'. w2 + ^2^2 
On integrating with respect to uj we can write the structure factor Sh{k,t = 0) as 



(12) 



7 



^.(fc,t^0)./5.(fc,.)^.^ + ^ 



(13) 

Recognising that the scahng form of Sh{k,t = 0) '^ fc^^^^"**, we notice that at cannot in general be determined 
from Eq.(n3). This is because the second term on the right-hand-side of Eq.(0) dominates at smaU momenta k 
provided an > 0, indicating that ah is indeterminate to this order of calculation. 

We turn now to the critical exponents in p. The single loop self-energy Sp(fc,cij) is given as shown in FigJ^a by 



I]p(fc, cj = Q) = -f,'J^J ^Gp{k - q, -n)Sn{q, n)q 
Inserting the expressions for Gp{k — q,Lu ^ fl) and Sh{q, ^) we find 
Ep(fc,w = 0) 



(14) 



-M 



dq 

2^ 



dn 

2^ 



1 



iil + \k — q 

This gives, on performing the integral over internal frequency fi, 

dq q^ 



2q^ 



n^ + g22h 



l+2ah 



q 



Sp(fc,cj = 0) = -Ai' f 



27r gi+2"/> \k - ql'^" + g^ 



In order to discuss this further in the context of Zp, we need to make a statement about ah and Zh- We have already 
obtained z^ = 1 in the foregoing and will now quote our numerical result for ah, viz. ah = 0.5. For small k the 
self-energy can then be written as 



Sp(fc,Lj) 



dq 1 j^ f dq 1 

2n{q + q^p) ^^ J 2tt {q + q^") {q + q^^p) 



We see from the above that ^p{k, 0), the relaxation rate for p fluctuations, is negative and finite as fc ^ 0, and we need 
to add a positive constant, Eq, to the self-energy (Eq > |Ep(/c — > 0)|) for regulatory purposes. This divergence in the 
relaxation rate, needing regulation, is reflected in the divergence we have encountered in our numerical investigations 
below; we have there followed an analogous procedure by introducing a numerical regulator which replaces divergent 
values of the transfer term by suitably deflned cutoffs (s). The resulting constancy of Ep implies Zp « for the 
regulated equations and will be used in the following. 
The correlation function Sp{k,uj) is given by (FigJ|b) 



Sp{k,uj) = 



1 



dq 

2^ 



dO. 



ik-q)^Shik-q,LO-n)Sp{q,n) 



(w2 + fc22:p) J 27r 7 27r 

The above integral will now be evaluated in the limit g — > and since Q ~ q^'^ for Sh we can replace 

Sh{k - q,cj - ri) ~ Sh{k,uj). 
Then using the scaling relation Eq.(P) we have 

1 f dq Cp 



Sp{k,uj) 



k Cp . f dq 



P-k^Sh{k,u;) 



1 



27r gl+2"p 



Anap (w2 + fc2^p)(a;2-hA:2^'>) 



(15a) 
(15b) 
(15c) 



where the last step follows from introducing a low er cutoff fco in the momentum integration over q. 
Using Eq.(p|) we have after integrating Eq.(15c) over uj 



S'p(fc,t = 0)-fc"(^+2"'') 



„l-2ah 



k^'p^k^h -I- /c^p) 



(16) 



Finally using Zp w we have 



Zh 

ctp — cth-\- -z — 1 for large k (l^a) 

ckp = ah — 1 for small k (17b) 

Given our numerical result of a^ — 0.5, the above predicts a negative ap, at small k. This is consistent with, and 
validates our assumption of, a cutoff ko which arises naturally as the wavevector separating the region of ap < (no 



infrared divergence) and Op > (infrared divergence prevalent) in Eqs.(7b) and (10c) 



More importantly, this non-trivial result for Up indicates that should we see numerical evidence of a negative Up 
for small wavevectors, we will have verified the existence of an asymptotic hypersmoothing in our model equations, 
which has an important bearing on sandpile surfaces in the continuous avalanching regime. This is discussed further 
in our concluding section. 

B. Numerical Analysis 

We focus now on our numerical results for Case B. The coupled eq uations in this section and the following one were 
numerically integrated by using the method of finite differences |22(| . Our grids in time and space were kept as fine- 
grained as computational constraints allowed so that our grid size in space Ax was chosen to be in the range (0.1,0.5) 
whereas that in time was in the range Ai (0.001, 0.005). Thus the instabilities associated with the discretisation of 
nonlinear continuum equations were avoided and convergence was checked by keeping At small enough such that the 
quantities under investigation were independent of further discretisation. Our results were also checked for finite size 
effects. In the calculations of this section we chose Dh = Dp = 1.0 and fJ. — 1 and our results were averaged over 
several independent configurations. We have calculated the exponents a and /3 and the corresponding error bars using 
the linear least square fit so that — (1 + 2/3) and — (1 + 2a) are given by slopes of the fitted straight lines. 

On discretising the equations Eqs.(6) we found once again the divergences that were previously observed in [||. 
These divergences are in our view a direct representation of the infrared divergence mentioned above, and we follow 
here a parallel course in regulating these via an explicit regulator. In earlier work Iq], a regulator was introduced 
which replaced the function fip\7h by the following: 

T = +1 for np{Vh) > 1 

= Hp{Vh) for - 1 < pp(yh) < 1 
= -1 for np{Vh) < -1 

In addition in this paper, we have introduced noise reduction to the regulated equations which has led to a more 
accurate evaluation of all our critical exponents. 

The Fourier transform Sh{k,t = 0) (Fig.P) is consistent with a spatial roughening exponent at ~ 0.501 ±0.007 via 
our observation of 



Sh{k,t = 0)^k 



-2.03±0.014 



and the Fourier transform Sh{x — 0, lu) (Fig.[l(]|) is consistent with a temporal roughening exponent Ph ~ 0.465 ±0.008 
via our observation of 



Sh{x = Q,uj) ^ LO' 



-1.93±0.017 



Hence Zh ~ 1.07, and thus the exponent a;^ ~ (Eq.(H)), indicating that the p, vertex does not renormalise. 
Using ah ^ 0.5 in Eq.(|T2) we can write the structure factor Sh{k,uj) as 



Sh{k,uj) 



l + f72(l + f72) 



1 + n^ 1 



r2fc2 mp 



(18) 



where J7 — uj/Tok. We find from the above that the expected form of Sh{k,uj = 0) in the limit of small wavevectors 
to be 

5;,(fc,W = 0)-fc-3 (19) 

Realising that our computed a/i < 1, we obtain from Eq.(|l2|) the prediction 

Sh{k^O,Lu)r^uj-^ (20) 



The full structure factor Shik,uj) has been calculated at two different k points and Fig.|Tl| displays our results to 
Eq.(|l8l). The sohd and the dashed line in the Fig.|ll| are the plots of the Eq.® for k = 0.1 and k = 0.2 with Tq = 
0.4 and 0.5 respectively. The spatial structure factor Sh{k,ijj — 0) shows a power-law behaviour (Fig.|l^) given by 

in qualitative accord with Eq.(ll9), and the temporal structure factor Sh{k = Ojw) shows a power-law behaviour 
(Fig.|l^) given by 



Sh{k ^{),uj) ^ Lo 



-1.91±.017 



in accord with Eq.(P0|) 

Given our values of au — 0.5 and Zh — 1, Eqs.(17a) and (17t) predict a crossover in ap from at large k to -0.5 as 



fc — > 0. The single Fourier transform Sp{k,t ~ 0) (Fig.|l^) shows a crossover behaviour from 

5p(fc,t = 0)~fc-2-i2±o.oi7 

for large wavevectors to 

Sp{k,t — 0) ~ constant 

as fc ^ 0. In Fig.|l^ we find a crossover from 0.56 at large k to -0.5 as fc ^ 0, which shows the same trend as 
the prediction above. Note however that the simulations also manifest in addition to the above the normal diffusive 
behaviour represented by ap = 0.56 at large wavevectors. The single Fourier transform in time Sp{x = 0,w) (FigJlJ) 
shows a power-law behaviour 



Sp{x — 0,Ll>) ^ LO' 



-1.81±0.017 



While the range of wavevectors in Fig. 15 over which crossover in Sp{k,t — 0) is observed was restricted by our 
computational constraints, the form of the crossover appears conclusive. Checks (with fewer averages) over larger 
system sizes revealed the same trend; additionally our theoretical calculations support the observed crossover via Eqs. 
(17). 

C. Homing in on the physics: a discussion of smoothing 

We focus in this section on the physics of the equations and our results. In the regime of continuous avalanching 
in sandpiles, the major dynamical mechanism is that of mobile grains p present in avalanches flowing into voids 
in the h landscape as well as the converse process of unstable clusters (a surfeit of Vh above some critical value) 
becoming destabilised and adding to the avalanches. Our results for the critical exponents in h indicate no further 
spatial smoothing beyond the diffusive; however, those in the species p indicate a crossover from purely diffusive to 
an asymptotic hypersmooth behaviour. Our claim for continuous avalanching is as follows: the flowing grains play 
the major dynamical role as all exchange between h and p takes place only in the presence of p. These flowing grains 
therefore distribute themselves over the surface filling in voids in proportion both to their local density as well as to 
the depth of the local voids; it is this distribution process that leads in the end to a strongly smoothed profile in p. 
Additionally, since in the regime of continuous avalanches, the effective interface is defined by the profile of the flowing 
grains, it is this profile that will be measured experimentally for, say, a rotating cylinder with high velocity of rotation. 

IV. ANOMALOUS SMOOTHING: THE CASE OF TILT AND BOUNDARY-LAYER EXCHANGE 

(CASE C) 

The last case we discuss in this paper involves a more complex coupling between the the stuck grains h and the 
flowing grains p as follows 

dhix, t) _ jj^^2^^^ t)-T + f]{x, t) (21a) 



dt 

DpV''p{x,t)+T (21b) 



dp{x,t) _ ^2 
dt ~ " 



T{h, p) = -i'(V/i)_ - Ap(V/i)+ (21c) 
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with ri{x,t) representing white noise as usual. 
Here, 

z+ = z for z > 

= otherwise (22a) 

2_ = z for z < 

= otherwise (22b) 

This equation was also presented in earlier work J5| in the context of the surface dynamics of an evolving sandpile. 
The two terms in the transfer term T represent two different physical effects which we will discuss in turn. The first 
term represents the effect of tilt, in that it models the transfer of particles from the boundary layer at the 'stuck' 
interface to the flowing species whenever the local slope is steeper than some threshold (in this case zero, so that 
negative slopes are penalised). The second term is restorative in its effect, in that in the presence of 'dips' in the 
interface (regions where the slope is shallower, i.e. more positive than the zero threshold used in these equations), the 
flowing grains have a chance to resettle on the surface and replenish the boundary layer g . We notice that because 
one of the terms in T is independent of p we are no longer restricted to a coupling which exists only in the presence 
of flowing grains: i.e. this model is applicable to intermittent avalanches when p may or may not always exist on the 
surface. In the following we examine the effect of this interaction on the profiles of h and p respectively. 

The complexity of the transfer term with its discontinuous functions precludes any attempts to solve this model 
along the lines of the earlier ones. We make some remarks here, however, on the likely critical behaviour of this model. 

We observe that the transfer term 

T=-Ap(V/i)+-i^(V/i)_ 

can be thought of as a formal infinite series by invoking a suitable representation for the Heaviside step functions in 
Eq.(23). We are then led to consider the following more general structure for the transfer term T, 

oo oc 

T = -Xp{\7h) - v{Vh) - Y^ i/„(V/i)"+i ^pY^ A„(V/i)"+i (23) 

n— 1 n=l 

Note however that this is not a very well-defined expansion because the coefficients in the infinite series could well be 
very large, if not infinite. However, given this disclaimer, we can still make the following comments in the spirit of 
self-consistency i.e. subject to numerical verification. 

If Ap(V/i) were the only nonlinearity, as in Case B, we would have z?i = I. Using h ~ x°'^ and p ~ x"'', we see 
Ap(V/i) is a more relevant nonlinearity than z/i(V/i)^, the leading nonlinear term in the expansion of (V/i)_, and is 
likely to be the controlling nonlinearity for the extreme long wavelength behaviour. Fig.^ shows that the A vertex 
never renormalises in the presence of the the KPZ term z/i(V/i)^, so that Zh is always fixed at unity. However, the KPZ 
vertex corresponding to z/i(V/i)^ has distinct behaviour in different wavevector ranges. In the range where the vertex 
renormalises, we cannot say much about the behaviour of ah] however, in the range where it does not renormalise, we 
might imagine that normal KPZ hyperscaling at + z^ = 2 would be restored. This, with Zh = 1, would give at — ^. 

If Zh — 1, we can write the scaling relation Sh{k,uj = 0) for the double Fourier transform at zero frequency as 

5,,(fc,C^-0)~fc-2-2a, 

which, in the regime where the KPZ hyperscaling holds, should look like Sh{k,uj = 0) ~ k~'^. 

We now try to obtain additional insights into the behaviour of these equations using the Hartree-Fock approximation. 
The spirit of Hartree-Fock is to replace nonlinear terms by linear ones with coefficients that are generally determined 
self-consistently. To undertake that here, we note that the step functions (Eq.ESh give rise to nonlinearities and hence 
the simplest thing to do is to replace them by an expectation value (the argument of the step function is a random 
variable and hence this is an acceptable approximation). We represent this expectation value by a number c with 
< c < 1. The equations of motion thus read 

dh 

— = DhV^h - X'pVh - ly'Vh + r]h{x, t) (24a) 

^ = DpV'^p + X'pVh + v'Vh (24b) 

ot 

with A' — cX and v' — (I — c)i^ and are identical to the ones studied by Bouch aud et a l. H. We expect at least in 



some regime of Eqs.(21) to reproduce the mean-field results appropriate to Eqs.(24a,24b) 
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1. Results for the single Fourier transforms 

The single Fourier transforms Sh{k,t — 0) (Fig.^ and Sh^x — 0,a;) (Fig.^ show power-law behaviour corre- 
sponding to 



Sh{k,t = 0) r^k 



Shix = 0,ll!) '^ uj' 



-2.56±0.060 



-1.68±0.011 



which implies that the roughness and the growth exponents are given by respectively a/,. — 0.78 ± 0.030 and Ph — 
0.34 ± 0.005. This suggests Zh = au/ Ph ~ 2 contradicting the prediction oi Zh = I by perturbative methods and 
suggesting that the mean-field approach outlined in the above might be more appropriate. We discuss this further in 
what follows. 

However, the small k limit of Sh{k,t — 0) indicates a downward curvature and thus a deviation from the linear 
behaviour at higher k (Fig.|l7|). This curvature, which had also been observed in previous work |5|, indicates a smaller 
roughness exponent ah there, i.e. an asymptotic smoothing. 

2. Results for the double Fourier transforms 

The double Fourier transforms Sh{k,uj — 0) (Fig.nOh and Sh{k = 0,lu) (Fig.EO) show power-law behaviour corre- 
sponding to 

5,,(fc = 0,c^)~^-i-80±0-™7 

Sh{k,uj = 0) ^ ^-4.54±o.08i Jqj. lojge wavevectors 
~ constant for small wavevectors 

The double Fourier transform Sh{k — 0, uj) shows the usual uj^^ behaviour that we have seen before in Eqs.(|^) and 
( po| ) which we have already discussed earlier. 

The structure factor Sh{k,U! ~ 0) signals a dramatic behaviour of the roughening exponent ah, which crosses over 
from 

• A value of 1.3 indicating anomalously large roughening at intermediate wavevectors, to 

• A value of about —1 for small wavevectors indicating asymptotic hypersmoothing. 

The anomalous roughening ah = 1 seen here is consistent with that observed via the single Fourier transform 
(Fig.|l^) and suggests, via the perturbative arguments given previously, that Zh = 1. However, if we assume Zh = 2 
according to the results of the single Fourier transforms given above, this would lead to an ah of about 0.8, in 
agreement with the values obtained both via single Fourier transforms in the present paper, and in [p|. In either case, 
our values of ah (either 1.3 or 0.8) suggest anomalous roughening of the interface at moderately large wavevectors. 

The anomalous smoothing obtained here (a^ '^ —1 if Zh '^ 1, and ah '^ —1.5 in the event that Zh is taken to be 2) 
is also consistent with the downward curvature in the single Fourier transform Sh{k, t — 0), as both imply a negative 
ah', we mention also that the wavevector regime where this smoothing is manifested is almost identical in both Figs.O 
and|l|. 

Since we expect that the anomalous smoothing results from a failure of the expansion of the st ep functions along 
the lines of Eq.(p3), this underlines our expectation that the mean- field solution of Eqs.(^4a],^4b|) would capture at 



least some of the flavour of this regime. We have therefore solved the mean-field equations (Eqs. |24a ,24b) numerically, 
and from Fig.E2and Fig.p3 we find that there is a crossover in Sh{k, t — 0)(Fig.p3) from a diffusive behaviour {zh — 2) 
at high wavevectors to a smoothing behaviour at low wavevectors. 

This behaviour is reflected in our results for Case C. At low frequencies the region of anom alous sm oothing can 



be understood by comparison with the corresponding region in the mean-fleld equations Eqs.(24a,24b) which also 
manifest this. At large fc, Sh{k,t = 0) and Sh{k,uj = 0) indicate anomalous roughening with ah ^ Zh ~ 1 which 
is consistent with the infrared divergence discussed in the previous section. However, as in Case A, Sh{x = 0,a;) is 
dominated by the diffusive z^ = 2 arising from the presence of S{uj — v'k) in the mean-fleld solution of Case C. This 
behaviour is corroborated by an evaluation of the full structure factor S{ki,uj) (Fig.pi|) which shows a distinct peak 
at an uii given by w,; — v'ki] this is reminiscent of the Lorentzian obtained in Case A (Fig.Q^). In fact, to leading 
order, Sh{k,uj) can be fitted to a Lorentzian; however, as we reduce the relative strength of v{Vh)- with respect 
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to Ap(V/i)+ the Lorentzian peaks disappear, and we begin to see a 'shoulder' reminiscent, as it should be, of the 
behaviour observed in Case B (Fig.|ll|) . This suggests that the present model is an integrated version of the earlier 
two, reducing to their behaviour in different wavevector regimes; we speculate therefore that there are two dynamical 
exponents (z^ = 1 and Zh = 2) in the problem. 

V. DISCUSSION AND CONCLUSIONS 

We have presented in the above a discussion of three models of sandpiles, all of which manifest asymptotic smoothing: 
Cases A and C manifest this in the species h of stuck grains, while Case B manifests this in the species p of flowing 
grains. We reiterate that the fundamental physical reason for this is the following: Cases A and C both contain 
couplings which are independent of the density p of flowing grains, and are thus applicable for instance to the 
dynamical regime of intermittent avalanching in sandpiles, when grains occasionally but not always flow across the 
'bare' surface. In Case B, by contrast, the equations are coupled only when there is continuous avalanching, i.e. in 
the presence of a finite density p of flowing grains. 

The analysis of Case A is straightforward, and was undertaken really only to explain features of the more complex 
Case C; that of Case B shows satisfactory agreement between perturbative analysis and simulations. Anomalies 
persist however when such a comparison is made in Case C, because the discontinuous nature of the transfer term 
makes it analytically intractable. These are removed when the analysis includes a mean-fleld solution which is able 
to reproduce the asymptotic smoothing observed. 

We suggest therefore an experiment where the critical roughening exponents of a sandpile surface are measured in 

1. a rapidly rotated cylinder, in which the time between avalanches is much less than the avalanche duration. 
Our results predict that for small system sizes we will see only diffusive smoothing, but that for large enough 
systems, we will see extremely smooth surfaces. 

2. a slowly rotated cylinder where the time between avalanches is much more than the avalanche duration. In this 
regime, the results of Case C make a fascinating prediction: anomalously large spatial roughening for moderate 
system sizes crossing over to an anomalously large spatial smoothing for large systems. 

Finally we make some speculations in this context concerning natural phenomena. The qualitative behaviour of 
blown sand dunes ||2^ is in accord with the results of Case B, because sand moves swiftly and virtually continuously 
across their surface in the presence of wind. By contrast, on the surface of a glacier, we might expect the sluggish 
motion of boulders to result in intermittent flow across the surface, making the results of Case C more applicable to 
this situation. It would be interesting to see if the predictions of anomalous roughening at moderate, and anomalous 
smoothing at large, length scales is applicable here. 
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APPENDIX: 



In this appendix we discuss some of the technical points related with the double Fourier transform. We have found that the 
crossover that we have seen in the Eq.(5) would not have been observed had we been using the single Fourier transforms 
Sh(k, t = G) and Sh{x — 0, uj) for numerical purposes. We illustrate this by writing explicitly the expressions for the relevant 
quantities: 



Sh(fc,i = 0)~fc-" 

Sh{x = 0,Lli) ^ UJ~ 
Sh{x = 0,tj) ~ L0~ 



for Lo small 
for LO large 



(Ala) 
(Alb) 
(Ale) 



The examination of Sh{k,t — 0) (FigJ4|) on its own yields no indication of the crossover to the smoothing fixed point; 
although there is a crossover in the Sh{x — 0,lu) graph (Fig.ga) from uj~^''' to cj~^, the analysis below shows that both 
regimes reflect diffusive behaviour, so that the smoothing fixed point {ah = 0,Ph = 0, Zh = 1) is entirely suppressed. 
The single Fourier transform Sh(x,UL)) is defined by 



Sh{x = 0,uj) 



dk 



Sh{k,uj) = 



dk 



Dhk''^ 



{uj ~ck)^ + Dlk^^h 



27r ' ' ' / __ 2-n- Dhk^h 
In the limit uj ^ ck the term in the square brackets behaves like a (5-function and thus 

" dk 



Sh{x = Q,L0)^ 



=-— ; 5(lO — Ck) ~ 

27r Dhk^'h io^h 



(A2) 



This is the origin of ballistic behaviour in the fiow term and is responsible for two anomalies. 

1. Firs tlv. w e notice from the above that the 5-function causes Sh{x = 0, ui) to behave like aj~^''. Comparing with 
Eq.( Alb) this leads to Zh ~ 2. However a simple-minded application of Eqs.(|3|) would have led to the uirong conclusion 
of (3h ~ 0.5. Even if the correct scaling relation Eq.(A2) were employed, the ballistic nature of the fiow term picks 
out, misleadingly, the high frequency (diffusive) dynamical exponent in the low frequency regime of Sh{x = 0,a;) 
(Eq.(Alb)). The low wavevector, low-frequency smoothing behaviour is thus entirely suppressed. 

2. Secondly, spurious oscillations are observed (Fig.Blb)) in the graph for Sh{x = 0,uj) as a function of grid size. A 
consideration of the form of the structure factor Sh{x = 0,^) makes it clear the crossover from small to large 
ui should not involve any imaginary quantities, and therefore strictly speaking we should not see any oscillatory 
behaviour in the structure factor in this limit. However the full form of the structure factor Sh{x, ui) for finite x does 
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contain imaginary portions, which are responsible for the oscillations. The characteristic length and time scales in 
our problem are given by 

to = Dh/c xo = Dh/c 

Whenever grid sizes in time or space are comparable to these characteristic lengths, the profile fiuctuates across 
these intervals, which is then aggravated by the shock fronts associated with the flow term. This results in oscillatory 
behaviour arising from the non-zero intervals in x associated with the sampling of the profile to generate the Fourier 
transform, Sh{x = 0,tj), which introduce a flavour of Sh{x,Lo) for finite x. These become increasingly violent as 
c increases because of the increased fluctuations associated with the ballistic flow term over the grids. In order to 
avoid these oscillations, one should choose grid sizes Aa; and At in such a way that they are always less than the 
characteristic scales in the problem, i.e., 

Ax < a;o and At < to- 

In view of the above, it is necessary to use the double Fourier transform to obtain an unambiguous picture of the structure 
factor and to pick out the asymptotic smoothing although this strategy might on flrst appearance seem to be a compu- 
tational overkill. The overwhelming advantage is that, by scanning the structure factor as a function of frequency uj for a 
fixed k, one immediately sets two frequency scales ck and D^k^, thus making it possible to pick up the relevance of these 
scales in Shikju). 

We also mention that our discussion is equally applicable to the Kardar-Parisi-Zhang (KPZ) equation pM with the addition 
of a flow term. Here too, the use of the double Fourier transform reveals the presence of the 'smoothing' flxed point due 
to the flow term. 

FIG. 1. The correlation function Sh{ki,ijj) against uj for three different wavevectors fci = 0.02(0), k2 = 0.08(+) and 
fcs = 0.12(D) with parameters c — 2.0, Du = 1.0 and A^ — 1.0. The positions of the peaks are given by ui\ — 0.04, lj2 ~ 0.16 
and ujo, = 0.24 as expected from Eq.(H). 



FIG. 2. The double Fourier transform, S{k,uj — 0) obtained from Eq.(Uaj) (Case A) for the h-h correlation function showing 
the crossover from high to low k. The different markers in the figure correspond to different grid sizes Aa; to sample distinct 
regions of k space; thus the markers A, x, and □ correspond to decreasing grid sizes and increasing wavevector ranges. The 
parameters used in the calculation are c = Dh = A^ — 1.0 and the characteristic wavevector is ko = tt- — 1.0. The dashed line 
is a plot of Sh(k,ijj = 0) vs A; for Case A with appropriate parameters, to serve as a guide to the eye. 

FIG. 3. The double Fourier transform, S{k = 0,tj) vs lo obtained from Eq. (|4a|) (Case A) for the h-h correlation function. 
The different markers in the figure correspond to different grid sizes At to sample distinct regions of tj space; thus the markers 
A, X, and □ correspond to decreasing grid sizes and increasing frequency ranges. The solid line is a plot of Sh{k = 0,^;) vs oj 
for Case A with appropriate parameters, to serve as a guide to the eye. The parameters are c — Dh = A'^ = 1.0. 

FIG. 4. Log-log plot of the single Fourier transform Sh{k,t = 0) vs fc obtained from Eq.(|4a|) (Case A) with parameters 
c = Dh = A^ — 1.0. The best fitted line shown in the figure is given by a slope of —1 — 2q;, — —1.90 ± 0.016. The characteristic 
wavevector fco is given by ko = tt- = 1.0 

FIG. 5. (a)Log-log plot of the single Fourier transform Sh{x = 0,uj) vs uj obtained from Eq.(Wa|) showing a slow crossover. 
Lines 1 and 2 in the figure are the best fits in the low and high ui regions with slopes —1 — 2f3h = —1.87 ± 0.003 and 
-1 - 2l3h = -1.525 ± 0.006 respectively. 

(b) Log- log plot of the single Fourier transform Sh{x = 0, lj) vs uj for two different values of c; c = 10 and c = 5 for data sets 1 
and 2 respectively. Note the increase in oscillation for increasing values of c. The other parameters are Dh ~ A^ = 1.0. 

FIG. 6. (a) The behaviour of < p{t) > as a function of time t. Here < p{t) > is the average over the sandpile surface of 100 
sample configurations. The grid size At = 0.005 and c — A = Dh = 1.0. 

(b) The root mean square width prms{t) — {< p^ > — < p >^)^ against time t over 100 sample configurations with parameters 
c = A^ = Dp = Dh = 1.0 

. (c) The variation of pmaxit) and pmin(t) with time t. pmaxit) and pmin{t) are respectively the maximum and minimum values 
of p for a given configuration of the sandpile at time t. Again, c = Dh = Dp — A^ — 1.0. 
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FIG. 7. One-loop diagrams for (a) Ep(fc,aj), the self-energy in p, (b) T,h{k,uj), the self-energy in h for the coupled equations 
of Case B (Eq.(6)). (c) The glossary for the diagrams shown in Fig.7(a,b) and Fig. 8. For example, the propagators for the h 
and p variables are represented by solid and dashed lines respectively, with a right arrow. Additionally there are diagrammatic 
definitions for the vertex and for the correlation functions for the h and p variables. 

FIG. 8. One-loop diagrams for (a) Sh{k,ij), the h — h correlation function, (b) Sp{k,uj), the p — p correlation function for 
the coupled equations of Case B (Eq.(8)). 



FIG. 9. Log-log plot of the single Fourier transform Sh{k,t = 0) vs fc obtained from Eqs.(6) ( Case B). The best fit has a 



slope of -1 - 2ah = -2.03 -I- ±0.014. Other parameters are p = Dh = Dp = /\i = Al = 1.0 



FIG. 10. Log-log plot of the single Fourier transform Sh{x — 0,1^) vs cu for Case B obtained from Eqs.(6). The best fit shown 
in the figure has a slope of -1 - Iolk = 1.93 ± 0.017. Again ^ = Dh = Dp = A| = Ap = 1.0. 

FIG. 11. The double Fourier transform Sh{ki,uj) vs lu { Case B ) calculated at two different wavevectors ki — O.l(O), 0.2(-|-). 
The curves ( solid(l) and dashed(2) lines ) shown in the figure are plots of Eq.([lg) with Fq = 0.4 and 0.5( for fei and k2 
respectively), to serve as a guide to the eye. Other parameters are p — 2, A^ = Ap = 0.1, Dh ~ Dp = 1.0. 

FIG. 12. Log-log plot of the double Fourier transform Sh{k,u! — 0) vs k { Case B ) obtained from Eqs.(6). The best fit has 
a slope of -(1 + 2ah + zn) = -3.40 ± 0.029. Again, p = 1.0, D^ = Dp = 1.0, A^ = A^ = 0.5. 

FIG. 13. Log-log plot of the double Fourier transform Sh(k — 0,ui) vs cu obtained from Eqs.(6) ( Case B ). The best fit 
displayed in the figure has a slope of —(1 + 2f3h H — —) ~ —1.91 ± 0.017. Other parameters are p — 1.0, Dh = Dp — 1.0 , 



FIG. 14. Log-log plot of the single Fourier transform Sp{x — 0,a;) vs uj obtained from Eqs.(6) ( Case B ). The best fit has 
a slope of -1 - 2/3p = -1.81 ± 0.017. Again, p = 1.0, Dh = Dp = 1.0, A| = A^ = 0.5. 

FIG. 15. Log- log plot of the single Fourier transform Sp(k,t — 0) vs k { Case B ) showing a crossover from a slope of 
-1 - 2q:p = at small k to -2.12 ± 0.017 at large k. Other parameters are p = 1.0, Dh= Dp = 1.0, A| = Ap = 0.5. 

FIG. 16. One-loop corrections to (a) the KPZ vertex, and (b) the A vertex for the coupled equations of Case C (Eqs.(21)). 

FIG. 17. Log-log plot of the single Fourier transform Sh{k,t — 0) vs fc for Case C obtained from Eqs.(21). 
The slope of the fitted line is given by —1 — 2ah = —2.56 ± 0.060 The parameters used in the simulation are 
u = m,\ = 1.0, Dh =Dp = 1.0, A| = 1.0. 

FIG. 18. Log-log plot of the single Fourier transform ^^(a: = 0,uj) vs u obtained from Eqs.(21) ( Case C ). The best fit has 
a slope of -1 - 2/3h = -1.68 ± 0.011 with parameters ;/ = 10, A = 1.0, Dh = Dp = 1.0, A| = 1.0. 

FIG. 19. Log-log plot of the double Fourier transform Sh(k,u) = 0) vs fc obtained from Eqs.(21) ( Case C ). The best fit 
for high wavevector has a slope of —(1 + 2ah + Zh) = —4.54 ± 0.081. As A: -^ we observe a crossover to slope of zero. Other 
parameters are Dh — Dp = 1.0, A^ — 1.0, v — IQ and A = 1.0. 
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FIG. 20. Log-log plot of the double Fourier transform Sh{k — 0,uj) vs oj obtained from Eqs.(21) ( Case C ). The best fitted 
line shown in the figure has a slope of —{l + 2/3h + — ) = —1.80 ±0.007. Other parameters are Dh = Dp = 1.0, A| = 1.0, i^ = 10 
and A = 1.0. 

FIG. 21. The double Fourier transform Shiki, uj) vs u obtained from Eqs.(21) ( Case C ) evaluated at three different 
wavevectors fci = 0.2(0), fc2 = 0.4(+) and k^ = 0.8(D) with parameters Dh = Dp = 1.0, A| = 1.0,;^ = 5 and A = 1.0. The 
peaks correspond to frequencies uj\ — 1.0, ci;2 ~ 2.0,0; = 4.0. 



FIG. 22. Log-log plot of the single Fourier transform Sh{k,t = 0) vs fc obtained from the mean-field Eqs.(24a 24b). 
The high k region is fitted with a line of slope —1 — 2ah ~ —2.05 ± 0.017. The low k region is fitted with a line of slope 
— 1 — 2ah = —0.93 ± 0.024. Note the crossover from ah ~ 0.5 at large k to zero at small k. Other parameters are v' — 10, 
A' = 2.0, Dh = Dp = 1.0, A| = 0.1. 



FIG. 23. Log-log plot of the single Fourier transform Sh{x — 0,uj) vs u for the mean-field Eqs.(24a,24h). The best fit has a 
slope of -1 - 2l3h = -1.94 ± 0.001 with parameters v' = 10, A' = 2.0, Dh = Dp = 1.0, A| = 0.1. 



17 



1e+08 



1e+07 - 



1e+06 - 



5 

^ 

5 



100000 - 



10000 - 



1000 - 



100 




0.001 



0.01 



0.1 
w 



10 



1e+10 



o 

II 

CO 



1 - 



1e-10 



Fig.2 




0.01 



0.1 



10 



100 



1000 



1e+10 



Fig.3 



5 

o" 

II 

CO 



1 - 



1e-10 




0.001 



0.01 



0.1 



10 
w 



100 



1000 



10000 



1 00000 



o 

II 

-1— • 

















Fig.4 




12 






















1 1 




1 




1 1 1 1 1 


















Data o 








o 

N 










Best fit — - 




11 


— 


% 












— 


10 


- 




~^^« 

X 

X 


•-^o 








- 


9 








\ 


"^^^^ 


^-^ 






8 














--o^ 

^~^-^<^ 
""^x 




7 


















6 


- 












V 


- 
































^.o 


















♦>^^ 




5 














o^4 




4 


- 












- 


Q 






1 1 




1 




1 1 1 1 1 





-5.5 



-5 -4.5 



-3.5 



ln(k) 



-2.5 



■1.5 -1 



5 

II 




Fig.5(b) 



5 

II 




ln(w) 



2.0e-06 



Fig.6(a) 



1.0e-06 - 



o'O.Oe+00 

V 



■1.0e-06 - 



-2.0e-06 




O.Oe+00 



2.0e+05 



4.0e+05 



Time (t) 



6.0e+05 



8.0e+05 



1 .Oe+06 



Fig.6(b) 



■g 

w 

E 
DC 




100 



1000 



10000 
Time(t) 



100000 



1e+06 



X 

CO 

E 



- 



"E 



■1 - 



1 || 



* i| 



Fig.6(c) 




« 
II 

,1' h 
II 

ii 



f I! 
' II 



I II 



II 



' !l 



ii I 



|!|i 

III, i,v 



I 

' ' ' 
,1 II 

i I I 



Al 



I' 
if 



!. 



1' 
II 
II 



'I 

I 'I 
|| I 

III 

I III 

i(||i 

ii 
\ 



^ 



* 
I' 

< I' 
I ii ii 

ii n 
ill I ' 

I 
11 (« 

< 
I 



III 

'i I'li'i/ 



ill I ' 
II 



i!iil 






ii»' 

III 

III 
i III 
I II II 



I I ii 
( I 



1 1, ii. 



III! 
U I I 



li' I "I I 



I l| 
Till 
l| 
l| 
l| 
'|l 
'I 



I I' 
I I 



,) 1 



- i 



Ij 



i' I II ,11 1 
•i 1(1 "ill 



1 1 



iiii 



n I I i 

* II I » I 

I II I II I 

I h ii II I 

I II ii II I 

I II i| II I 

I II I'll' I 

I" '!! ! 



II II' 1 

III Ml I 
III III! I 
M 11] 
III II I 
111 II 

III II 

I" " 
I' ' " 
I"" 
|l I II 
(11 
I II ' 
I I II 
I < 
I 



O.Oe+00 



2.0e+05 



4.0e+05 



Time(t) 



6.0e+05 



8.0e+05 



1 .Oe+06 



^p(k,(o) 



— fc,co 




^,co 



k— q, (i) — Q. 



(a) 



q,Q 



-O-:- 



.^,^ 



\(kM 



— fc,co 




V ^,co 



k— q, (i) — Q. 



(b) 



-^-^-"- Gp(A:,ca) 


- -KX- - Sp(A:,CO) 


^ ~ ^ — r r/r m^ 




'J/jIa, wj 


q,Qy^ 




^x — 


|lp(V/?j vertex 






(c) 




Fig. 7 






k— q, Od — Q. 



k— q, (i) — Q. 



(a) 




k— q, Od — Q. 



k— q, (i) — Q. 



(b) 
Fig. 8 



Fig.9 



o 

II 
-1— • 




Fig.10 



5 

II 




ln(w) 



1000 



CO 



100 - 



10 - 



1 - 



0.1 - 



Fig.l1 



line 1 — 

k1 o 

line 2 --- 

k2 + 




0.01 



0.01 



0.1 



w 



10 



Fig.l2 



5 




5 

II 




14 



o 



Fig.l4 



~i 

Data o 
Best fit — 



12 - 



■^ 



>. 



5 

o" 

II 

o 



10 - 



8 - 



6 - 



4 - 



2 - 



.o 



o~~ 



o 








-4 



ln(w) 







Fig.l5 



o 






h 




+ 




Dressed KPZ vertex 



Bare KPZ vertex 



(a) 




.-/ 



h 




.-/ 



h 



/ P 



+ 




Dressed 'kp(\/ h) vertex 
Bare Xp(\/h) vertex 



(b) 



Fig. 16 



26 



24 - 



22 - 



20 - 



Fig.17 

— I 



Data o 
Best fit — 



o 

II 
-1— • 



CO 



"I 



18 - 



^ 16 - 



14 - 



12 - 



10 - 



o 



o 



o 



o 




8 



-8 



-6 
ln(k) 



-4 



14 



13 - 



12 - 



Fig.18 

1 



♦ 



^. o 



Data o 
Best fit — 



^ 



11 - 



-^ 



^-^ 



5 

II 



10 - 



9 - 



8 - 



7 - 



6 - 



-6 






*oX 



Oh 







-4 
ln(w) 



-2 



Fig.l9 



5 




5 

II 



20 



18 - 



16 - 



14 - 



12 - 



10 - 



8 



Fig.20 




-4 



ln(w) 







10000 



1000 - 



<> 



d 100 It 



10 r 



[3 



0.01 



o 



+ 



Fig.2l 



o 



<> 



o ^o 



o 



o 






+ 



"^ + + 4- + 



D Q ° □ □°°Qn°n^n°#^' 



0.1 



k1 o 
k2 + 

k3 □ 




10 



w 



16 



14 - 



12 - 



Fig.22 



T 



Data o 
line 2 — 
line 1 --- 



o 

II 
-1— • 



o 



10 - 



8 - 



6 - 



4 - 



>^-^. 



■^- 



;-^-.4 



~o^ 



o-<^o 



^ %^ 




-4 



ln(k) 







Fig.23 



5 

II 




